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ABSTRACT 

, Context. The coupling between time-dependent, multidimensional MHD numerical codes and radiative line emission is of utmost 

• , importance in the studies of the interplay between dynamical and radiative processes in many astrophysical environments, with 

i-C ' particular interest for problems involving radiative shocks. There is a widespread consensus that line emitting knots observed in 

Herbig-Haro jets can be interpreted as radiative shocks. Velocity perturbations at the jet base steepen into shocks to emit the observed 
spectra. To derive the observable characteristics of the emitted spectra, such as line intensity ratios, one has to study physical processes 
<^ , that involve the solution of the MHD equations coupled with radiative cooling in non-equilibrium conditions. 

^ ' Aims. In this paper we address two different aspects relevant to the time-dependent calculations of the line intensity ratios of forbidden 

^ ' transitions, resulting from the excitation by planar, time-dependent radiative shocks traveling in a stratified medium. The first one 

concerns the impact of the radiation and ionization processes included in the cooling model, and the second one the effects of the 
numerical grid resolution. 

Methods. Dealing with both dynamical and radiative processes in the same numerical scheme means to treat phenomena characterized 
^ ' by different time and length scales. This may be especially arduous and computationally expensive when discontinuities are involved, 

1/^ I such as in the case of shocks. Adaptive Mesh Refinement (AMR) methods have been introduced in order to alleviate these difficulties. 

fSJ , In this paper we apply the AMR methodology to the treatment of radiating shocks and show how this method is able to vastly reduce 

^SJ ■ the integration time. 

Results. The technique is applied to the knots of the HH 30 jet to obtain the observed line intensity ratios and derive the physical 
parameters, such as density, temperature and ionization fraction. We consider the impact of two different cooling functions and 
different grid resolutions on the results. 

Conclusions. We conclude that the use of different cooling routines has effects on results whose weight depends upon the line ratio 
considered. Moreover, we find the minimum numerical resolution of the simulation grid behind the shock to achieve convergence in 
the results. This is crucial for the forthcoming 2D calculations of radiative shocks. 
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1. Introduction that differ by orders of magnitude, and this becomes a major 

problem when tackling the time-dependent evolution of such ra- 

Supersonic flows are ubiquitous in the Universe: expanding su- diative shocks. For example, in Herbig-Haro jets, one needs to 

pernova remnants, stellar winds, AGN and Herbig-Haro jets are treat in the same scheme spatial scales wefl below ~ 10'^ cm, 

some examples, and shocks are abundant and prominent in these to reproduce the time-dependent post-shock temperature varia- 

flows. Shocks located in extragalactic environments, such as tions correctly, and scalesJ^lO^^ c m to stu dy the behavior of the 

AGN jets, can be considered adiabatic, since the cooling time electron density (Massaglia et al. I2005al Paper I). Under these 

for thermal emission typically exceeds by far the source lifetime, conditions and with these requirements, employing an Adaptive 

On the other hand, shocks that form in galactic objects such as Mesh Refinement (AMR) technique becomes almost mandatory, 

supernova remnants and Herbig-Haro jets must be treated in- This technique can provide adequate spatial and temporal reso- 

cluding radiative effects. lution by dynamically adapting the grid to the moving regions of 

Radiative shocks have been studied in ste ady-sta te condi- interest, giving a tremendous saving in computational time and 

tions by several authors (e.g. Cox & Raymond [UMl Hartigan memory overhead with respect to the more traditional uniform 

et al. 119941 ). who derived the one-dimensional post-shock be- S^^'^ approach. 

havior of various physical parameters (temperature, ionization The important issue of numerical spatial resolution in time- 
fraction, electron density, etc.), as functions of the distance from dependent simulations of radiative shocks in 2D has been con- 



the shock front. From these studies it turns out that the phys- sidered by |Raga et al. (2007) employing an AMR technique as 



ical parameters vary behind the shock front on scale lengths well. They discussed the dependence of the morphological struc- 
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ture of the perturbations depending on the grid resolution. Their 
maximum grid resolution was with a cell size corresponding to 
1 .5 X 10'^ cm. They found that while the detailed structure of the 
shocks depends on the resolution, the emission line luminosities, 
integrated over the volume, are less dependent on cell size. 

A second crucial question is: does the dependence on the 
cooling function details have qualitative or quantitative ef- 
fects on the calculated distributions of the physical parameters? 
Te§ileanu et al. ( I2008I I have developed a detailed treatment in- 
cluding an ionization network of 29 species and compared, at 
constant numerical resolution, the resulting distribution of phys- 
ical parameters (such as temperature, density, etc.) with the one 
obtained employing the simplified cooling of Paper I, for a 2D 
cylindrical jet affected by perturbations that evolved in internal 
working surfaces. They found quantitative differences but quali- 
tatively the outcomes were very similar 

A more challenging problem than the determination of the 
integrated line emissions and distributions of physical parame- 
ters is the calculation of line intensity ratios between different 
species. In this case, the fractional abundances of the various 
species vary in very different ways with space behind the shocks 
(see Paper I). This variation is mainly governed by the temper- 
ature profile, which is typically a very steep function of space. 
Therefore we will examine the effects of both the details of the 
cooling function and of the grid resolution adopted. 

As discussed in Paper I and Massaglia et al. ( |2005bl ), obser- 
vations of some HH jets at distances of a few arcseconds from 
the source typically show that the behavior of temperature, ion- 
ization and density along the jet is incompatible with a freely 
cooling jet. It is now generally accepted that the line emission 
in HH jets is a result of the observation of several post-shock 
regions of high excitation with a filling factor of about 1%. In 
our calculations, we will refer, in particular, to the observational 
results described in Bacciotti et al. (1999*), where a specialized 
diagnostic technique has been applied to HST data of the HH 30 
jet. Hartigan & Morse (120071 1 also re-examined the HH 30 case 
using slitless spectroscopy performed with the Hubble Space 
Telescope Imaging Spectrograph. Their results are fully consis- 
tent with the findings in Bacciotti et al. ( I1999D . 

In this paper, we will follow the dynamical evolution of an 
initial perturbation as it steepens into a (radiative) shock travel- 
ing along the jet, and derive the post-shock physical parameters 
consistently. From these parameters, we construct the synthetic 
theoretical emission line ratios to be compared with observa- 
tions. Using this setup, we implicitly consider that the medium 
in the real jet returns close to the steady jet conditions between 
the propagating shocks. This was verified in parallel 2D simula- 
tions (Te§ileanu et al. 2009) , as the ID simulation, being done in 
the reference frame of the mean flow, lacks the steady flow of the 
jet. After the shock and the high-density post-shock zone pass, 
there follows an underdense zone, and then the medium returns 
to the stationary conditions of the constant flow. The distance 
after which this happens increases during the evolution of the 
shock from 0.5 x 10''*cm close to the jet base to about lO'^cm 
at 500AU from the base. This increasing distance becomes re- 
solved observationally at about 200AU from the base of the jet, 
and also at such distances interactions between different shock- 
waves become likely, these being possible explanations for the 
oscillations visible in the observational data after 200AU. The 
emission at each point of the jet is computed by simulating the 
propagation of a shock produced at the base of the jet up to that 
point. As discussed before, we will adopt two different treat- 
ments for the radiative losses, namely: i) the simplified cooling 
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Fig. 1. Snapshot of the temperature behavior vs distance for the 
uniform grid, for a propagating shock in a stratified medium. 



employed in Paper I and ii) the new cooling treatment described 
in Tegileanu et al. ( I2008I I. 

The plan of the paper is the following: In Section 2 we 
present the computation setup and the adopted techniques to 
model the problem; in Section 3 we apply the model to the case 
of HH 30 and discuss the results obtained with the two cooling 
functions and the effects of the grid resolution; the conclusions 
are drawn in Section 4. 



2. The model 

In this Section we give a general outline of the model and the 
form of the initial perturbation. More details on the setup can be 
found in Paper I, where we analyzed the case of the DGTau jet. 

As mentioned in the Introduction, in the preshock condi- 
tions assumed here (decreasing density departing from a few 
10'*cm"-', uniform temperature of 2000K) the post-shock tem- 
perature drops by about one order of magnitude over a scale 
length lower than ~ lO'^cm i.e. a much smaller distance than the 
post-shock evolution that goes on over distances of a few times 
10'^ cm. This has been a compelling reason for developing and 
employing Adaptive Mesh Refinement, as in Paper I, described 
in detail in a companion paper (Mignone et al. 2009 ). We have 
verified the validity of this approach by computing several of the 
models presented in the next section on uniform and adaptive 
meshes and found a tremendous gain in efficiency. Figs. ([T]i and 
^ show the results of one particular model computed on a uni- 
form grid employing 49152 zones and on an equivalent adaptive 
grid with 5 levels of refinement (base level resolution is 1536 
zones). On a Pentium IV processor with a 1.7GHz clock and 1 
GB of RAM, the uniform grid approach took about 6.04 x 10"* 
sec, while the AMR computation required only 238 sec, up to 
the same final integration time f = 5 x 10^*. 
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Fig. 2. The same as in Fig. [T]but for the AMR integration and 
five levels of refinement, for an equivalent resolution of 49152 
cells, with a zoom on the shock region. The histogram in the 
lower part of the plot (with the right vertical scale) represents 
the refinement level employed in each cell of the simulation. 



where Kisa constant derived from the single shock assumption, 
■y - 5/3 (the gas is monoatomic; molecules form at lower tem- 
peratures, while we are interested in the post-shock regions of 
high radiative emission), Uo is the mean flow velocity and the 
initial velocity perturbation is: 



u(x) 



uo[—(x — xi)^ + 2cr(x — xi)] if 2cr + xi > x > xi 
otherwise 



where mq is the perturbation amplitude, X[ is the initial coordi- 
nate of the perturbation and cr is its half-width. However, the 
exact shape of the perturbation is not crucial for the model. In 
our calculations we take xi - lO'^ cm and cr - 2x lO'"' cm. 

In carrying out the calculations we have set Uo{x) - 0, i.e. 
the reference frame is that of the mean flow (that is, the steady 
jet flow on which the perturbation is set). To transforme the re- 
sults back to the laboratory frame we have set Uq = 100 km s"', 
which is the bulk velocity we assume for the HH 30 jet (Uo must 
not be mistaken for mq, the initial velocity perturbation ampli- 
tude). 

The boundary conditions assume free outflow at jc = and 
X = L. The extent of the computational domain, L, has been 
chosen sufliciently large to follow the shock evolution for f ~ 15 
yrs and to avoid spurious interactions with the boundaries. For 
this reason, we adopt L = 4.5 x 10'^ cm. 

The assumption of a monoatomic gas is justified in the case 
of HH30 and many other YSO jets - however, there are cases 
(the so-called "molecular jets") when this assumption is not valid 
anymore, and molecular cooling must be taken into account. 



2.1. Model equations and initial perturbation 

We restrict our attention to one-dimensional MHD planar flow, 
which implies that we are following the shock evolution along 
the jet axis. The fluid is described in terms of its density p, ve- 
locity M, thermal pressure p and (transverse) magnetic field By. 
These quantities obey the standard MHD equations for conserva- 
tion of mass, momentum, magnetic field and total energy, in the 
presence of radiative cooling represented by the energy loss term 
J1{T, X) (energy per unit volume per unit time), which depends 
on the temperature T and on the ionization state of the plasma, 
described by the vector of ionization fractions X. A discussion 
on the cooling function follows in the next subsection. 

A nonuniform preshock density that decreases away from the 
star is a reasonable assumption when dealing with an expanding 
jet, as suggested by observations. We note that shocks propa- 
gating into a stratified medium tend to increase their amplitude 
when they find a decreasing preshock density. We consider the 
following form for the preshock density: 



po{x) = po 



x'^+xP 



(1) 



where x is the spatial coordinate along the jet axis, xq sets the 
initial steepness of the density function (this affects the shock 
evolution even at larger distances) and I < p <2. 

As in Paper I, the form of the initial perturbation has been 
taken in such a way as to keep the hydrodynamic Riemann in- 
variant J =constant. In this way a single forward shock is 
formed. In this scheme, the density perturbation is: 



r ■ 



1 



(u-Uo) +P(,' 



(2) 



2.2. Radiative cooling 

In the current work we use two different approaches for the com- 
putation of the radiative cooling loss term X.(T, X). In a first sim- 
plified model (see Paper I), X consists of the ionization fraction 
of H only, whereas a second, more accurate treatment is summa- 
rized below. 

In the detailed cooling treatment, described in Tegileanu et 
al. (120081) , 28 additional evolutionary equations are solved for 
the non-equilibrium ionization fractions X (of H, He, C, N, O, 
Ne and S). The loss term accounts for energy lost in lines as well 
as in the ionization and recombination processes. Line emissions 
include contributions from transitions of the 29 ion species. The 
metal abundances are the ones adopted by Bacciotti et al. ( 19991 1, 
in particular A^/// = 1.1 x 10"^, 0/H = 6 x 10"^ S/H =1.6 x 
10"^ and the abundance of C is 10% of the solar one. 

To carry out comparisons between the observed and com- 
puted line intensity ratios, we have determined the populations 
of the atomic levels relevant for [SII], [Nil] and [OI] emission, 
solving the excitation - de-excitation equilibrium equations for 
five energy levels, according to ,Osterbrock & Ferland (2006)] 
The line emissions reported in Bacciotti et al. ( 1 19991 1 for the 
HH 30 jet are: [SII]]/16716 + ^6731, [NII]^6548 + ^6583 and 
[OI]/16300, and the corresponding intensity ratios [OI]/[NII] and 
[SII]/[OI]. 



3. The case of the HH 30 jet 

To illustrate the astrophysical application of the above method- 
ology, we will apply it to the numerical simulation of some of 
the physical properties of the HH 30 jet. 

We have integrated the magneto-fluid equations employing 
the AMR method described previously with eight levels of re- 
finement, using the PLUTO code (Mignone et al. 120071 ). As in 
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Paper I, in order to obtain values directly comparable with ob- 
servations, we have space averaged all post-shock quantities at 
each evolutionary time point: 



(2(0) = 



JQ(x, t)e{[SII](x,t)}dx 
J e{[SII](x,t)]dx 



(3) 



where Q(x, t) is a physical quantity such as electron density or 
ionization fraction or temperature. This procedure is applied be- 
cause by processing the observational data, we only have access 
to the physical parameters of the line-emitting regions, so this 
weighting is implicit. Differently from Paper I, where we av- 
eraged the line intensity ratios using the total line emissivity of 
the corresponding emitters as weighting functions, in the present 
paper we find the unweighted average of the line emissivity and 
calculate the ratio of the average lines afterward. This procedure 
does not lead to results that differ qualitatively from the previ- 
ous ones, but appears closer to the actual observation process. 
To calculate the line emissivities, the atomic transition rates and 
collision strengths adopted are the ones described in Tegileanu et 
al. (2008 ). The collision strengths have been interpolated using a 
Lagrangian scheme to account for their temperature dependence. 

The main free parameters of the model are: the particle den- 
sity (at X = 0) riQ, xo, the velocity perturbation amplitude mq 
and initial transverse magnetic field Bq. Taking advantage of 
the peculiar efficiency of the AMR method, we have widely ex- 
plored this parameter space to find good agreement with the ob- 
served Une ratios. We have set the (uniform) preshock tempera- 
ture Tq = 1, 000 K and the initial ionization fraction /; = 0.1%, 
i.e. due to the metal contribution only. 

Both radiative cooling treatments were employed in the sim- 
ulations, searching for the best agreement of the model results 
with the observations for each one. A discussion of the differ- 
ences in the results follows. 

3.1. Simplified cooling 

The simplicity of this cooling model allows for rapid simula- 
tions, and thus the very efficient exploration of the parameter 
space. The root mean square deviations of the simulated line ra- 
tios with respect to the observational ones were computed (a few 
examples are shown in Table [TJ in order to quantitatively esti- 
mate how close each simulation is to observations. The simula- 
tions presented in the table have base densities of 2 x lO'^cm"^ 
(n2e4) or 5 x 10'*cm"-' (n5e4), transverse magnetic field 300//G 
(b300) or 500yuG (b500), and amplitude of the perturbation in 
velocity 50 (v50), 55 (v50) and 70km s"' (v70). 



Simulation 


[OI]/[NII] 


[SII]/[OI] 


n5e4b500v70 


0.440 


0.305 


n2e4b500v55 


0.161 


0.103 


n5e4b500v55 


0.121 


0.143 


n2e4b500v50 


0.068 


0.110 


n2e4b300v50 


0.063 


0.108 



Table 1. Root mean square deviations of line ratios with simpli- 
fied cooling. 



A combination of the parameters that yielded a reasonably 
good agreement with the observed line intensity ratios, employ- 
ing the simplified cooling model, was hq - 2 x 10'* cm"^. 
Bo - 300//G, Xo - 0."1 and mo - 50 km s In Paper I 
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Fig.3. Intensity ratios of [SII]]^6716 + /16731, [NII]^6548 + 
/16583 and [OI]/16300 lines vs distance. Symbols represent the 
data from Bacciotti et al. (1999) and curves are the computed 
model curves (using the simplified cooling model): [OI]/[NII] 
(bullets, dot-dashed line) and [SII]]/[OI] (circles, solid line). The 
model parameters are no = 2x 10"* cm"-'. Bo = 300//G, xq = 0."1 
and Mo = 50 km s"', with Uo = 100 km s"' as the jet bulk veloc- 
ity. 



for the DG Tau jet we set p - 2, in the present case a better 
agreement is obtained setting p - 1, i.e. for a less steep de- 
crease of the preshock density with distance, i.e. a slower expan- 
sion rate for the jet as it propagates into the ambient medium. 
This can be explained by the higher degree of collimation of 
the HH 30 jet. A comparison with Paper I shows also that mq is 
lower in this case with respect to the DG Tau parameters, i.e. 
Mo = 50 km s"' instead of mo = 70 km s"'; this is reflected in 
the value of the ionization fraction, considerably lower for HH 
30. Moreover, we note a higher value of Bo that must be as- 
sumed in the present case with respect to DG Tau. The results 
are shown in Fig. [3] where we plot the observed line intensity 
ratios of [OI]/[NII] (bullets) and [SII]/[OI] (circles) (as reported 
by Bacciotti et al. 1999 ) and, superimposed, the corresponding 
model curves (([OI]>/([NII]> dot-dashed line and <[SII])/<[OI]) 
solid line). Since the model proposed is extremely simple, we 
mainly aim, as shown in Fig. [3] to obtain values of the inten- 
sity ratios and trends with distance that are reasonably close to 
the ones observed, and we cannot account for the knots of emis- 
sion visible along the HH 30 jet, resulting in ample spatial os- 
cillations of the line intensity ratios, especially in the [OI]/[NII] 
ratio. 

Note also that a different choice of the metal abundances 
can shift the model lines of Fig. |3] by some percent; e.g. 
ILodders (2003)| reports N/H = 0.802 x 10"^ 0/H = 5.81 x lO""* 
and S/H= 1.83 x 10^1 

In Fig. |4] we plot the physical parameters that are de- 
rived with a special diagnostic technique from observations by 
Bacciotti et al. (1 19991 ) (symbols) compared to the outcome of the 
present model (lines). We note that the electron density derived 
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Fig. 4. Average physical parameters of the jet derived from ob- 
servations following Bacciotti et al. ( 1 19991 1 (symbols) and from 
the model with simplified cooling (curves). Temperature (dia- 
monds, solid line), electron density (circles, dashed) and ioniza- 
tion fraction (bullets, dot-dashed). 



from the shock model (dashed line) exceeds the one derived by 
Bacciotti et al. ( 119991 1 (circles) for distances below ~ 1". As dis- 
cussed in Paper I, however, when the electron density is above 
the critical density, the ratio of observed [SII] lines, from which 
rtg is derived, saturates, and the diagnostic in this condition, as 
reported in the aforementioned paper, yields only lower limits to 
the electron density. The discrepancy in the electron density af- 
fects the diagnostics of temperature (diamonds) that differs from 
the model one (solid line). 



3.2. Detailed cooling 

The previous cooling model might be an oversimplified approx- 
imation in the case of moderately powerful shocks such as the 
ones we are dealing with here. This leads to one of the motiva- 
tions of this paper: a comparative study of the differences be- 
tween and advantages of each of the two approaches to radiative 
cooling treatment. 

Employing the detailed cooling model, we have explored the 
parameter space for the same setup as in the previous case. The 
same method of estimating the root mean square deviation with 
respect to observations was employed. The best agreement with 
the observed line intensity ratios was obtained for the following 
parameters: no - 5x 10* cnT^, Bq = 500;uG, xq = 0."1, mq = 55 
km s ' and p - I. The resulting line ratios are plotted in Fig. 
|5]and the evolution of the physical parameters in Fig. |6l In this 
scheme, the full ionization state of the plasma is computed at 
each step, allowing for non-equilibrium states. This of course 
also has a serious impact on the simulation speed, because of the 
supplementary equations in the system. 

The spatial profile of the physical parameters is similar to the 
one obtained in the best agreement case employing the simpli- 



Fig.5. Intensity ratios of [SII]]^6716 + /16731, [NII]^6548 + 
/16583 and [OI]/16300 lines vs distance. Symbols represent the 
data from Bacciotti et al. (1999) and curves are the computed 
model curves (with detailed cooling treatment): [OI]/[NII] (bul- 
lets, dot-dashed line) and [SII]]/[OI] (circles, solid line). The 
model parameters are no = 5x 10"* cm"-', Bq = 500;uG, xq = 0."1 
and Mo = 55 km s"', with Uo = 95 km s"' as the jet bulk veloc- 
ity. The dashed, lighter lines are the results obtained with the 
simplified cooling model for the same parameter set. 
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Fig. 6. Average physical parameters of the jet derived from ob- 
servations following Bacciotti et al. (1999) (symbols) and from 
the model with detailed cooling (curves). Temperature (dia- 
monds, solid line), electron density (circles, dashed) and total 
ionization fraction (bullets, dot-dashed). 
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fied cooling. The same considerations of the underestimation of 
the electron density close to the source apply. 

By applying the set of simulation parameters resulting from 
the simulations with detailed cooling to the setup employing 
simplified cooling, important differences appear. Due to the 
higher overall density of the plasma and the much simpler chem- 
istry in the simplified cooling case, the electron density remains 
at values above the ones derived from observations by a factor 
of 3 along the evolution of the shock. This is also reflected in the 
line ratios that depart from the observed ones by a factor of 2 for 
the [01] /[Nil] and 20% for the [SI I]/ [01] ratios, triggered by 
the higher maximum temperatures attained. 

3.3. On the importance of numerical resolution 

It is important to understand the resolution needed for these 
simulations, before moving to 2D configurations (Te§ileanu et 
al. 120091 1. A rough estimation of the resolution needed in or- 
der to resolve the post-shock zone can be done as follows: for 
a maximum peak temperature of about lO^K and densities of 
< 10^ cm"^, the internal energy density of the gas is of the or- 
der of lO^^ergcm"^. According to previous computations (see 
Te§ileanu et al. 120081 ), the maximum energy loss rate in radia- 
tive cooling processes in these conditions is ~ 10"'^ erg cm^^s ' . 
Allowing for a maximum relative change in internal energy of 
1 % in a timestep, the size of the timestep is restricted to about 
lO'^s, which combined with the jet velocity of a few hundreds 
kms~' leads to a space scale of ~ 10^'cm. Even higher resolu- 
tions, of a few times 10'"cm, are needed to accurately resolve the 
ionization/recombination processes (with timescales of < lO-'s). 
For the setups used in this work, this spatial resolution is attained 
with 8 levels of refinement. 

We have studied the effect of increasing resolution on the 
line ratios and found it to be important up to a saturation value 
above which the line ratios do not change significantly. For this, 
we have computed the root mean square deviations between the 
simulated values of the line ratios with an increasing number of 
refinement levels. The results are presented in Table |2l where 
the first column shows the number of refinement levels in the 
simulations between which the RMS deviations were computed. 
Each refinement level doubles the resolution. Three configura- 
tions are presented: one with base density N=2 x 10'*cm"'', trans- 
verse magnetic field 300juG, perturbation amplitude in velocity 
50km s"'; the second one has N=5 x lO'^cm"^^, B=100yuG, ve- 
locity perturbation 40km s"' ; the third one has N=5 x lO'^cm"-', 
B=500/iG, velocity perturbation 55kms The convergence of 
the results can be inferred from these data. 



Sim. 


n2e4b300v50 


n5e4bl00v40 


n5e4b500v55 


ref.levs. 


Rl 


R2 


Rl 


R2 


Rl 


R2 


base ~ 2 


0.1378 


0.0319 


0.1240 


0.0477 


0.1150 


0.0477 


2-4 


0.0512 


0.0171 


0.0883 


0.0534 


0.0884 


0.0499 


4-6 


0.0141 


0.0816 


0.0425 


0.0418 


0.0589 


0.0509 


6-8 


0.0110 


0.0141 


0.0088 


0.0145 


0.0256 


0.0359 


8-10 


0.0004 


0.0008 


0.0017 


0.0023 


0.0033 


0.0084 



Table 2. Root mean square deviations of line ratios (Rl is 
[OI]/[NII], R2 is [SII]/[OI]) obtained with simulations of in- 
creasing resolution. 



Presented in Fig.|2]are four cases of increasing resolution for 
the first configuration in Table|2](n2e4b300v50): one employing 
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Fig. 7. Effect of the simulation resolution on the line ratios of 
[SII]]^6716 + ^6731, [NII].16548 + /i6583 and [OI]^6300 lines, 
in the case of simplified cooling. The setup is identical to the one 
in Fig. |3] with the four line sets representing results from simu- 
lations with 10 levels of refinement (long-dashed lines), 6 levels 
of refinement (dash-dotted lines), 3 levels of refinement (dashed 
lines) and without refinement at all (dotted lines), respectively. 



only the base grid (1 536 zones), a second one with three lev- 
els of refinement (equivalent resolution 12288 cells), one using 
six levels of refinement (equivalent resolution 98 304 cells), and 
one with 10 levels of refinement (equivalent resolution 1 572 864 
cells). It can be seen that the line ratios converge to the results of 
Fig. [3] obtained in a simulation with 8 levels of refinement. The 
simulations presented in Fig.Quse the simplified cooling model, 
but similar effects were obtained when using the detailed cooling 
model. In Fig. [T] the line ratios obtained employing 8 levels of 
refinement almost perfectly overlapped the ones obtained with 
10 levels of refinement, so they were not shown. 

The physical spatial resolution in a simulation with 8 levels 
of refinement in the present setup is about 1.2 x 10 '"cm, and 
should be enough to adequately resolve the post-shock zone, 
with its rapid variation of physical parameters. The resolution 
for 6 levels of refinement is ~5xlO'°cm, it drops to 4 x 10"cm 
for three levels of refinement and to 3.25 x 10'~cm at the base 
grid resolution. 

The line ratios are very sensitive to resolution variations, 
making it important to reach an adequate (high) resolution for 
reliable results. This is a critical aspect for future 2D simula- 
tions, that will have to employ AMR to reach these extremely 
high resolutions with the available computational power. 

In order to see how the resolution requirements depend on 
the parameters of the shock, in Fig. |8]we present the evolution 
of the temperature at the same age for three Shockwaves with 
different amplitudes of the initial perturbation (decreasing from 
left to right). The x axis does not represent the actual position 
of the shocks but only a distance scale, the front of the shock- 
waves being plotted close to each other for convenience. We see 
that with the increase of the perturbation amplitude (and thus the 
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Fig. 8. Three simulations of Shockwaves with density N=2 x 
10'*cm"^, transverse magnetic fields 300/iG, and perturba- 
tion amplitudes in velocity of (from left to right) 70, 50 
and 40km s"', respectively, at the same evolutionary point. 
Logarithmic plots of temperature. 
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Fig. 9. Three simulations of Shockwaves with transverse mag- 
netic fields 300//G, perturbation amplitude in velocity of 
50km s"\ and base densities of (from left to right) 8, 5, 
and 2 x lO^cm"-', respectively, at the same evolutionary point. 
Logarithmic plots of temperature. 



shock strength), the evolution after the shock steepens. In the 
same way, in Fig.|9]three Shockwaves with different initial den- 
sities are plotted, with decreasing density from left to right. The 
evolution of the temperature after the shock steepens when in- 
creasing the density. Typically, as the dependence steepens, the 
resolution needed increases. All the cited plots were made for 
the same average velocity of the underlying jets. When decreas- 
ing this velocity, the spatial resolution needed for an accurate 
analysis of the post-shock zone increases. 

3.4. Comparison and discussion 

The parameter sets for the best agreement obtained in the previ- 
ous paragraphs were not the only ones that well interpret the ob- 
servations, but other possible "successful" combinations do not 
depart from these values of the parameters. As an example, set- 
ting Uo = 95 km s"' for the bulk jet velocity, we obtained model 



curves very similar to those of Figs.|5]and|6]with the same p - I, 
XQ = 0."1 but with «o = 4 X lO'' cm"^ Bq = 300//G, and mq = 58 
km s ' . The same considerations apply for the simplified cooling 
simulations. 

The jet bulk speeds Uq were chosen in such a way that the 
total shock speed in the rest reference frame would be approxi- 
mately equal to 150km s"' in all cases - this being an approxi- 
mate average knot speed derived from observations of Doppler 
shifts and proper motions (Hartigan & Morse |2007] |. 

Comparing the results obtained with the two approaches to 
radiative cooling, identical or similar values of the model param- 
eters can be inferred (Table |3]l. Parameters such as velocity am- 
plitude of the perturbation (determining the shock strength) and 
transverse magnetic field (resulting in shock compression) are 
similar in the two cases. The total density of the jet is however 
one of the most difficult parameters to estimate from observa- 
tions, and the different results in our two cases come from the 
different evolution of the total radiative losses and ionization in 
the post-shock zone. 



Parameter 


Simplified 


Detailed 


no (cm-') 


2x10^ 


SxlO-* 


fio(pG) 


300 


500 


Mo(km s"') 


50 


55 


f/o(kms-') 


100 


100 




0.1 


0.1 


7'o(K) 


1000 


1000 



Table 3. Comparison of the model parameters that gave the best 
agreement with the observations of HH30. 



The simplified radiative cooling treatment appear to be 
suited for a rapid exploration of the parameter space, giving re- 
sults (parameters of the model) close to the ones obtained with 
more sophisticated cooling. 

4. Conclusions and summary 

We have demonstrated in previous works that numerical inte- 
gration of time-dependent radiating (shocked) flows can sub- 
stantially benefit from the employment of adaptive mesh refine- 
ment methods, and the present work supports this conclusion. 
By computing the cooling time scale according to the algorithm 
described in Mignone et al. (2009), one could take full advantage 
of the time adaptation process. The use of this method allows us 
to tackle simulations of radiating shocks in 2D as they propagate 
along a cylindrical supersonic jet. 

Our results applied to the DG Tau jet (Paper I) as well as 
the HH 30 jet (present paper) revealed a good agreement be- 
tween the observed and calculated line intensity ratios of differ- 
ent transitions and for different elements (Figs.[3]and|5]l, and also 
of derived physical parameters, including the ionization fraction 
(Figs.|4]and|g. 

An accurate treatment of the radiative losses allows us to 
more realistically compute the full ionization state of the plasma, 
and thus accurate line emission for each ion, in the approxima- 
tion of a five-level atom. This, however, requires significant ad- 
ditional computing power that is not always available. An alter- 
native is the use of the simplified cooling model, following only 
the ionization state of hydrogen, that proved to give good indi- 
cations on the perturbation amplitude and transverse magnetic 
field of our model. This simpler but much faster strategy may be 
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used for a preliminary exploration of the parameter space, being 
followed by simulations with detailed cooling only in the areas 
of interest. Both cooling models were discussed and compared 
in Te§ileanu et al. (2008 ), and in the present work we extended 
the comparison to line ratio diagnostics and application to a real 
HH object. 

Grid resolution proved to be a critical parameter for the line 
ratios estimations from MHD simulations. The physical dimen- 
sions of the computational zones have to not exceed about 10" 
cm in order to accurately capture the evolution of the physi- 
cal parameters and emission properties in the post-shock zones. 
This is a crucial aspect for future 2D simulations that will have 
to fulfill these requirements. This result can be generalized for 
the typical conditions encountered in shocks traveling through 
protostellar jets, as the variability in these conditions could not 
change the numbers by an order of magnitude. Weaker shocks 
(shocks evolving from smaller initial perturbations) or lighter 
shocks (lower initial density) will be less stringent on the resolu- 
tion requirements, while stronger and/or denser shocks will add 
a factor to the present results. 

We can summarize the general picture as follows: i) when 
one is interested in the distribution in space of the jet physi- 
cal parameters such as density, temperature, etc., the details of 
the assumed cooling function matter very little (Tegileanu et al. 
12008 ) provided the numerical resolution suffices to minimize nu- 
merical dissipation effects; ii) when one considers the detailed 
shock structure, numerical resolution is more important, but less 
so for the behavior of the integrated emission line luminosities 
( Raga et al. (2007)} ; iii) when one discusses the line intensity 
ratios, a finer resolution must be achieved, while the details of 
the cooling function have effects that differ according to the line 
considered, in the present case from about 10% up to a factor of 
two. 
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